% Script to load HER data and plot two types of splines
% 1) annual 365 days knots 2) seasonal 45 days knots on magnetic data at a
% station. This is done on request by Stefan (See mail on Oct 29, 2014)
% The idea is that there can be long departures of seasonal variations in magnetic data 
% from the annual splines.

close all
clear

% load data

ncfname = '/Users/manojnair/data/obs_mag_data/WDC_Intermagnet/update_2013/absolute/HRN_1995_2012.nc';
start_index = 0; % start from the beginning
ndata = 9468000; % read all data
time_array_min = (datenum(1995,1,1,0,0:9468000-1,30));
[x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj] = read_geomag_netcdf(ncfname, start_index, ndata, 0);

 f_data = sqrt(x_data.^2 + y_data.^2 + z_data.^2);
     declination_est = nan(size(x_data));
        
        L = isnan(x_data) | isnan(y_data);
        
        declination_est(~L) = 180/pi * atan( y_data(~L)./ x_data(~L) );

% read abs x values



L = isnan(f_data);
data_array = f_data(~L);
time_array = time_array_min(~L);

b1 = min(time_array):365:max(time_array)+10;
b2 = min(time_array):45:max(time_array)+10;
    
    sp = spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
    v = ppval(time_array_min,sp);
    sp2 = spline(b2,data_array(1:60:end)'/spline(b2,eye(length(b2)),time_array(1:60:end)'));
    v2 = ppval(time_array_min,sp2);
figure(1);
set(gcf,'Position', [560   369   787   579]);
    
subplot(311);
%plot(time_array_min,f_data);
set(gca,'FontSize',16);
hold on;
plot(time_array_min(~L),v(~L),'r','LineWidth',2);
plot(time_array_min(~L),v2(~L),'g','LineWidth',2);
ylabel('F');
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
datetick('x','keeplimits');
h=legend('data','annual fit','seasonal fit');
set(h,'FontSize',8)
%% declination

L = isnan(declination_est);
data_array = declination_est(~L);
time_array = time_array_min(~L);

b1 = min(time_array):365:max(time_array)+10;
b2 = min(time_array):45:max(time_array)+10;

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
    v=ppval(time_array_min,sp);
    sp2 = spline(b2,data_array(1:60:end)'/spline(b2,eye(length(b2)),time_array(1:60:end)'));
    v2 = ppval(time_array_min,sp2);

subplot(312);
%plot(time_array_min,declination_est);
set(gca,'FontSize',16);
hold on;
plot(time_array_min(~L),v(~L),'r','LineWidth',2);
plot(time_array_min(~L),v2(~L),'g','LineWidth',2)
ylabel('Declination');
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
datetick('x','keeplimits');



clear  v data_array time_array L;

%% read abs z vlues

L = isnan(x_data);
data_array = x_data(~L);
time_array = time_array_min(~L);
%fit spline


b1 = min(time_array):365:max(time_array)+10;
b2 = min(time_array):45:max(time_array)+10;

    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
    v=ppval(time_array_min,sp);
    sp2 = spline(b2,data_array(1:60:end)'/spline(b2,eye(length(b2)),time_array(1:60:end)'));
    v2 = ppval(time_array_min,sp2);


subplot(313);
%plot(time_array_min,x_data);
set(gca,'FontSize',16);
hold on;
plot(time_array_min(~L),v(~L),'r','LineWidth',2);
plot(time_array_min(~L),v2(~L),'g','LineWidth',2)

ylabel('X');
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
datetick('x','keeplimits');
% set(gcf, 'PaperPositionMode', 'auto');
% saveas(gcf,'HRN_1','epsc');


%% Make plot of a subset
subplot(311);

axis([datenum(2000,1,1) datenum(2003,1,1) 5.405e4 5.425e4]);
h=legend('annual fit','seasonal fit');
set(h,'FontSize',8)

subplot(312);
axis([datenum(2000,1,1) datenum(2003,1,1) 2.75 3.75]);
subplot(313);
axis([datenum(2000,1,1) datenum(2003,1,1) 7900 8050]);

set(gcf, 'PaperPositionMode', 'auto');
saveas(gcf,'HRN_2','epsc');

%% plot without originl data
subplot(311);
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
subplot(312);
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
subplot(313);
axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
subplot(311)
set(gcf, 'PaperPositionMode', 'auto');
saveas(gcf,'HRN_3','epsc');

